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Abstract 

We define a sparse hermitian lattice Dirac matrix, H, coupling 2n + 1 Dirac 
fermions. When 2n fermions are integrated out the induced action for the last 
fermion is a rational approximation to the hermitian overlap Dirac operator. 
We provide rigorous bounds on the condition number of H and compare them 
to bounds for the higher dimensional Dirac operator of domain wall fermions. 
Our main conclusion is that overlap fermions should be taken seriously as 
a practical alternative to domain wall fermions in the context of numerical 
QCD. 
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I. INTRODUCTION 



A major embarrassment of lattice field theory in the context of QCD has disappeared: we 
now have a way to preserve chiral symmetry on the lattice , and the all important the- 
oretical understanding of "soft" physics consequences of chiral symmetry at the Lagrangian 
level can be taken over from the continuum to the lattice. Naturally, one is eager to exploit 
this development in numerical QCD, and at the moment there are two ways that have been 
explored: One is the so called domain wall fermion approach P]|3|§], and the other is based 
on implementing ||] a rational approximation to the sign function of the overlap Dirac oper- 
ator [[J . The relation between the two is elucidated in [|7| . One conclusion from |7| could be 
that with computers immensely more powerful than the ones we have at present it wouldn't 
matter which approach one uses for numerical QCD. In practice, the two methods are quite 
different, and it is important to assess their relative strengths and weaknesses. 

To compare methods one needs to separate the quenched case from the dynamical case: 
In a dynamical hybrid Monte Carlo simulation one needs to invert the Dirac operator often, 
at each step of the so called "trajectory" . The inversions do not generate propagators used 
in computing physical observables, except when the trajectory is completed and the change 
in gauge fields accepted. In a quenched simulation one extracts physics results from each 
set of fermion propagators. Thus, the comparison of overall fractional costs for inversions 
works out differently in the quenched and in the dynamical case. 

It seems that the overlap approach has some advantage in the quenched case: it has been 
possible to do simulations at very low quark masses ||, lower than what was achieved with 
domain wall fermions ||. To be sure, a systematic and complete comparison has not been 
carried out, so this is more an impression than a hard fact. The reason for the difference is 
that with the overlap it has been possible to better handle the cases where the Hermitian 
Wilson Dirac operator, Hyy, has eigenvalues very close to zero. Numerical investigation 
has shown that there is a finite density of eigenvalues very close to zero at typical lattice 
couplings used in lattice simulations . While this is a problem for both domain walls |7j 
and overlap, in the overlap case one can exploit the simple local structure of Hy/ and project 
out the troublesome states ]TI|. This is expensive, but needs to be done only once per gauge 
field configuration and will help in the calculation of all propagators at this gauge field. 
Thus, the extra expense is amortized in the quenched case. A similar projection method in 
the case of domain wall fermions might make both methods equal. Until very recently there 
existed no practical implementation of a projection technique in the domain wall context, 
but this is changing as we write ||12|| . 

At the dynamical level it always was felt that domain wall fermions were superior because 
the action was quite standard in form and one needed only one conjugate gradient inversion, 
rather than the two nested ones required by the rational overlap Here again, this 



is more an impression than a hard fact because one ought to take into account condition 
numbers, matrix sizes and the possibility to use projectors. But, superficially at least, it is 
hard to ignore the advantage of domain wall fermions with a straightforward action and a 
relatively well tested algorithm. 

The latter distinction between the implementation of domain wall fermions and rational 
overlap fermions could be eliminated however by undoing what the overlap does relative 
to domain wall fermions. One reintroduces extra fields which interact quadratically by a 
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sparse hermitian matrix H. The main requirement of H is the following: There exists one 
massless field, such that integrating out all the other fields produces for ijj a specific 
effective action S B &. S c r = —ipH^ip and the approximate hermitian overlap Dirac operator 
is given by if™ = 5(75 + e n (H w )) where the function e n (x) is a numerically accurate rational 
approximation to the sign function sign (a;) for x in the spectrum of Hw- One can easily add 
an explicit Dirac mass term for the field ip. 

H has to fulfill additional requirements: It has to have a condition number that is not 
significantly worse than that of the five dimensional Dirac operator used for domain wall 
fermions, D. In particular, the condition number should not diverge as a function of the 
number of extra fields, 2n. The cost of acting with H should grow no more than linearly 
with n, so long as the truncation of the sign function converges as e~ cn . Furthermore, H 
should depend structurally only on Hw as one entity. A dependence on Hw as one object 
ensures that if a better version of Hw is eventually found, it would be trivial to change 
H, replacing the old Hw with the new one. Also, the implementation of the projection 
technique simply involves a replacement of Hw as a whole. In particular, taking derivatives 
with respect to the gauge fields would be simple, making the computation of the "force" in 
hybrid Monte Carlo relatively easy. Although the action for domain wall fermions also has a 
simple dependence on the gauge fields, the dependence on the length of the extra direction 
is controlled by the complicated transfer matrix Tw 0. The dependence on n is simpler in 
the overlap it comes in only through e n (Hw), and Hw is a sparse matrix, unlike 

Tw- 

Previous proposals for H had some numerical entries that were large and grew with 
n. This implied that the norm of H, \\H\\, grew and made the condition number k(H) = 
1 1 if 1 1 1 |if _1 | I grow too. Our choice for matrix norms is quite standard: | \X\ | = J A max (Xt X) , 
with A max (A^A) being the largest eigenvalue of X'X. The main objective of this paper is to 
present a new version of H which satisfies the above requirements and also has a reasonable 
condition number. 

In the next section we shall construct H. We proceed with the derivation of rigorous 
upper bounds on \\H\\ and on Rigorous bounds are nice, but there always is a 

question whether they are saturated. This will be discussed. We shall also ask the opposite 
question: what is the best one can hope for, regarding k(H) ? To answer that we shall derive 
a lower bound for This is all we need because, in practice, we expect the upper 

bound on \\H\\ to be typically almost saturated. We proceed to address the question how 
our results compare to what is known about the domain wall fermion operator D. Since not 
much seems to be known, we derive some exact results for D too, but also leave more work 
for the future. We then briefly discuss projection techniques for both overlap and domain 
wall fermions. 

We conclude the paper with a discussion of other advantages we envisage our proposal to 
have over domain wall fermions, beyond better chirality properties. We compare the number 
of Hw operations needed to perform the inversion of the fermionic operator. We find this 
number numerically for three cases: domain wall fermions, a method using the operator 
H proposed in this paper, and a method using a direct computation of the sign function 
represented by the same rational approximation as implemented by H. The comparison is 
carried out for the two dimensional Schwinger model. We hope to convince the reader that 
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the overlap alternatives merit serious numerical testing in the context of QCD, something 
we are not fully geared up to do efficiently by ourselves. 



II. CONSTRUCTION OF H 



The kernels of the quadratic fermionic actions we shall work with will be hermitian 
matrices in most cases. Let ip be the light Dirac field representing a light quark. We wish 
to end up with an action 



"75 



e n (H) 



w , 



(1) 



The bare quark mass /i is restricted by < 1 for physical reasons. 

The matrix Hw has the standard Wilson form, but can be easily replaced by a more 
elaborate construction. The function e n (x) is an approximation to sign(x): 



1 " 1_ 

n hi c 2 M w + 



with 



cos 6L 



7T 1 

sin6L 9 S = — (s ), s = 1, 2, n. 

2n T 



(2) 



(3) 



This rational approximation can be replaced by others: One can replace the c 2 s and s s 
quantities with other real numbers of either sign; to change the overall sign of a contribution 
at some s one simply switches the sign of Hw at that s. 

For each s we introduce a new field Xs- The \s fields are decoupled from each other and 
enter quadratically in the action: 



Xs 



c 2 s H w + 



H 



w 



Xs 



(4) 



In addition, they couple to the ip field by 



1 — /i 
2n 



J2(^Xs + Xs^)- 



(5) 



To get the right effective action for ip, there also is a quadratic term in i/j: 

^^75^- 



(6) 



To make the action local we should eliminate the inverse of Hw from the x-action. 
To this end we introduce n more Dirac fields, 4> s , also Grassmann, and change the x se lf 
interaction to: 

Xs{c 2 s H w )Xs + s s (Xs(ps + 4>sXs) - 4> s (H w )4>s 
Introduce the combined field ^ = (/0, xi, (fii, Xn, 4>n)- The action is 
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(7) 



with 



H = 





.fa 
V 2n 





V 2n 









Si 


—Hw 


J±=H 

V 2n 


















V 



1-M 
2n 















V 2n 



«2 











-H w 






2n 








c 2 n H w 










-Hw ) 



(8) 



Our new extended overlap model is based on the following identity: 



/ 



d^d^e 



ndet(c^ + ^) 

.s=l 



(9) 



The prefactor can be canceled by adding pseudofermions, which will be decoupled in the s 
index. 

At n = oo we can take a continuum limit, with 



de = —ds, o < e < -, 

2n ' 2' 



(10) 



The effective action then becomes 
1 +/i 



5 = 



2 + 



/i 



+ 



IT 



cos 2 9x(8)Hwx(0) 



Hw<p(0) + s\n9x{0)(t){0) + sin 



X 



The following identity (which holds as long as \\H W || is finite) expresses the essence of our 
construction: 



2 n 



-I 

7T JO 



d6 



2 r°° dt 



H w cos 2 6 + ^ 



/" 

7T JO 



£ ^W" _|_ t 

t Hw 



siga(Hw). 



(12) 



In terms of 9 the action in the extra dimension involves no derivatives. In a rough sense the 
6 variable corresponds to a "fifth" momentum (for a four dimensional physics application). 
The physical degrees of freedom, ip, couple only to an average field. These two features 
indicate that one might be able to use some sort of multi-grid or hierarchical techniques to 
reduce the 0(n) computational cost per inversion of H to a logarithmic dependence on n. 
Other integral representations of the sign function, or other approximation to the integrals 
above, will produce more variations on the same basic idea, but with possibly different 
computational properties. 
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The strategy for finding bounds on the spectrum of H is based on an exact formula for 
the determinant of H — z. This determinant is obtained by adding — z^f^f to the action and 
doing the gaussian integral. 



det(H - z) = I J] det [(c 2 s H w - z)(H w + z) + s 

ls=l 



det 



— ^— 75 + Z + ——f n {H W , Z) 



Here, 



f n (H w , z) 



1 

" ^ r;7/ 



(13) 



(14) 



Eigenvalues of H are roots of the equation det (if — z) = 0. All the roots come from roots 
of the last factor. (Roots of the factors in the product over s are canceled by poles in the 
last factor. So, the spectrum of H is determined by the last factor.) 
Introducing 



S n (a, b) 



we see that 



1 n 1 
Ly i 

n fr[ ac 2 s + bs 
1 



fn(Hw, z ) — S n (H w — z, 



z). 



(15) 



(16) 



S n (a, b) is a ratio of polynomials in a, b: 

2n 
2s- 



Q n (a,b) 



b n ~ s a s 



2n N 



- 1 , Q n (a,b) = J2[ 2s )b n - s a s . (17) 



s=l 



s=0 



Extracting (H w + z)~ n from numerator and denominator, we see that S(a, b) is also a ratio 
of polynomials in Hw- For real a and b we can write simple closed formulae: If ab > we 
have 



S n (a,b) = Slg "ffl tanh(ncu), oo = log f |1 + Jajb\/\\ - Ja/b\) ■ 
yab V / 



(18) 



If a& < 0, 



S n (a,6) = ^ffitan(ncu), < < J, = (1 + iJ-a/b)/{l - iJ-a/b). (19) 



These formulae make the n-dependence explicit. 
The product ab clearly plays a central role: 

ab = 1 + zg(H w , z), g(H w , z) = z — H w — 



Hw + z 



(20) 



Our bounds are based on the observation that for ab > we have an n-independent 
bound: 

1 



\f n (H w ,z)\ = \S n (a,b)\< 



ab 



(21) 



If we allowed ab < 0, we would have had little control over S n (a, b) because of the tangent 
function. 
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III. UPPER BOUND ON \\H\ 



For very large \z\ ab is large and positive, making z + ^-^f n (Hw,z) large in absolute 
value. It is clear that det(if — z) cannot vanish then. To find an upper bound for \\H\\ 
we look for the smallest \z\ values for which the above is still true, as we decrease \z\ from 
infinity. 

Theorem I: For any hermitian bounded Hy/ we have 



| [#11 < 

Proof of Theorem I: Pick a real number z, 

\z\ > \ \\H W \ 



Hw\\ 2 + 2. 



+ 2. 



(22) 



(23) 



We shall prove, by contradiction, that this implies det(if — z) ^ 0. All we need to show is 
that the hermitian operator ^-^75 + z + ^-^f n (H w , z) has no zero eigenvectors. We assume 
that there exists such a normalized eigenvector ip . ip obeys: 



1 + /i 



75^0 = - 



Z + 1 - - f n (H W , z) 



"00- 



(24) 



Taking norms we get 



1 + fi 



>4 



z + 



1-fX 



fn(Hw, z) 



^0 



Let h be an arbitrary eigenvalue of Hy/\ we know that 



|*| > Vh 2 + 2. 
Simple analysis of the function zg(h, z) implies 

zg(h, z) > 0. 

As a consequence: 

\fn(h,z)\<l. 

Hence, every eigenvalue A of z + f n (Hw, z) obeys 



(25) 



(26) 



(27) 



(28) 



IAI > \z\ 



1 — fi 



(29) 



(Recall that < 1.) By the variational principle, the right hand side of equation 
obeys the same inequality. Hence, 



1 + n 



> \z\ 



2 ~ 1 1 2 

which implies \z\ < 1 in contradiction to our initial assumption about z 



(30) 
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This establishes the upper bound of Theorem I. Note that it is /x and n independent. 
Corollary: Let Hy/ be the hermitian Wilson Dirac operator in even d dimensions, and H 
defined as above, with the standard replacement of the four dimensional 75. Let the mass 
parameter m in H\y be restricted by m > —2. Then 

\\H\ \ < \jl + (2d + m) 2 . (31) 

This is a direct consequence of the known upper bound on ||-£fV|| I 
We expect this bound to be quite close to optimal. 

This proves that the main difficulty faced by previous proposals to implement the overlap 
Dirac operator by adding extra fields has been completely eliminated. The problem with 
the previous proposals was that one could not control \\H\\. Moreover, our bound is quite 
stringent numerically, even smaller than a typical bound on a Wilson Dirac operator in 
2d+ 1 dimensions |14[]. This is a plausible comparison, because one could thing about s as 
indexing an extra dimension. At any rate, for d — 4 and m = —1.8 one would have 

||#||<6.4 (32) 

which is very reasonable. 

We ran some tests in two dimensions with m = —1 and found the bound of 3.3 to be 
typically almost saturated for £7(1) gauge configurations generated with a Wilson action 
at (3 = 2. Moreover, with a trivial gauge background one can explicitly check that one 
can get quite close numerically to the upper bound and so it is impossible to find a gauge 
background independent bound which is significantly better. 



IV. UPPER BOUND ON 1 1 1 1 1 

The immediate question we need to address now is what happens to the low eigenvalues 
of H. Having established that H behaves more or less as a usual fermionic lattice operator 
at high eigenvalues, numerical problems can only come from a large \\H~ ||. 

The basic strategy to get the bound is similar to the one used above. First take z = 
and n = 00. Assume that Hw has no zero eigenvalue. Then, ab — 1 and \\foo(Hw, 0)|| = 1. 



As a result, -^75 H — ^fooiHw, 0) cannot have zero eigenvalues if \i 7^ 0. We therefore look 
now for a neighborhood of z = where, if /i > and H\y has no zero mode, we have 



1 + ^ ^ 11 ^ l ~ ^ f lu \ 



(33) 



If this is true. 



^Is + z + ^UH^z) 



cannot have a zero mode. Indeed, if there were 



such a zero mode, ipo, we would have a contradiction, since then 



1 + (i 



^0^0 = % 



z + 



1 2 



f n (Hw, z) 



^0, 



(34) 



in violation of the variational principle for the maximal eigenvalue of 



z + 



w,z) 
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The range < fi < 1 covers all possible positive quark masses. The case — 1 < fi < 
describes quarks with a negative mass. This certainly is not without interest, but the analysis 
becomes more complicated and unnecessary for our purposes here. So, we simply restrict 
ourselves to the range < \i < 1. 

Let h be an eigenvalue of Hw- 

1 377 < \h\ < \\H W \\. (35) 



\\ H w II 

We define z n (h, //) as the smallest positive solution to the following equation: 

1 | 

— — = \z n (h,fi) H — f n (h, z n (h, (36) 

Since f n (—h, —z) = —f n (h,z), we have: 

z n (h,fi) = -z n (-h,/i). (37) 
h is restricted to ranges that are symmetric about zero. Therefore, the inequality 

1 -^>\z+ 1 -^f n (h,z)\ (38) 

can be guaranteed by restricting z to a range symmetric about zero: 

\z\ < u n (H w ,n). (39) 
Here, the numerical bound u n (H w ,fi) is given by 

u n (H w ,fi)= min z n (h,/j). (40) 

7^=T7<\h\<\\H w \\ 

\\H W | 

We have therefore proven the following bound: 

Theorem II: For < fi < 1 and arbitrary nonsingular hermitian Hw, we have: 

iij > u n (H w ,n). (41) 



\\H 

u n (H w ,fi) is determined by equations fl3"B"D and (|^) ■ 

z n (h, fj,) can be found numerically for the ranges of h and fi of interest in practice. To 
find a bound analytically is cumbersome because u n (H w , /i) will be controlled by either an h 
corresponding to a maximal (in absolute magnitude) eigenvalue of Hw or to a minimal one. 
Some analysis shows that there is a constant of order unity (dependent on the dimension d 
and the m parameter in Hw), C, and that if we restrict Hw by 

WHw'W^C, (42) 

the minimum over h we need in order to get u n (H w ,n) is attained at 

h = h min = (43) 
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Thus, the entire dependence of u n (Hw,^) on Hw comes in through the value of the lowest 
eigenvalue of H w . Gauge configurations for which ll-f/yj/ 1 !! < C are easily handled by any 
numerical method and the condition number in that case is not a source of concern. So, 
nothing is lost by assuming that the smallest eigenvalue of H w is smaller than C~ 2 ; this is 
anyhow the case for most gauge configurations one encounters in practical QCD simulations. 

In order to get some feeling for orders of magnitude we work out a simple analytical 
approximation for the bound on \z\ and represent the exact result for the bound as the 
product of this approximation times a correction factor. Some typical values for the cor- 
rection factor are presented in Table |. All the entries in Table | are for the case when the 
smallest eigenvalue of H w is smaller than C~ 2 . 

We now take ab > because this certainly is true at z = 0. As explained above, we can 
restrict our attention to positive z, without loosing generality. It is true now that 

\z + ^f n (h min , z)\\<z+ 1 ~ fi 1 (44) 

'1 + zg(h min , z) 



So, our bound will hold if we enforce the first inequality below: 



1 + u I — u 1 1 — it 

— — — > z + — ' 



+ zg(h min , z) 



2 



1 + * -= + 



h 



mm , 



1 — a 1 — a z 1 — a z . . 

" + ^ + V^"^ + 2T- (45) 

£j ii ''mm " z '"'min 

In the above series of approximations we assumed h m { n << 1. We end up with an approxi- 
mate range from which eigenvalues of H are excluded: 

\z\ < 2fih min . (46) 

Finally, we rewrite our exact result in the following form: 

iitj-iii > U n (H w ,fl) = C ra (ffy,^) — Ayr . (47) 
\\ H II \\ H w\\ 

For H-f^il > C the c n prefactor depends on H w just through h m ™ ■ In Table | we collected 
several values of c n (/i m in, A 4 ) to give a sense of the dependence on n, h min and \i. c n (Hw,(i) 
is calculated using equations (36), (40) and (47). 

V. CONDITION NUMBER OF H: WORST AND BEST CASE 

Combining our exact results we obtain an exact bound for condition numbers. 
Theorem III: 



<H) < v f tt \ (48) 
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TABLE I. The entries are numerically obtained exact values for the coefficient c n (h, /i) appear- 
ing in our rigorous bounds, h 2 always is the smallest eigenvalue of . 
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From Table | we see that decreasing n improves the condition number. This is expected 
because finite n effects induce chirality violation even at /i = 0, so act somewhat as an 
effective additional mass. We also see from Table | that at high enough n the entire correction 
needed to make the approximate bound rigorous is of order unity. The dependence of the 
correction factor on h and \i shows that the approximation is somewhat pessimistic. 

In summary, roughly, the main message is that the condition number of H cannot be 
worse than the product of the condition numbers for H\y and for the effective action gov- 
erning the light fermions, 

* (^75 + l -^e{H w )) ~ i. (49) 

It is natural to ask now: How pessimistic is the bound when viewed as an estimate for 
k(H) ? We shall answer this question in two ways: First, we shall find an upper bound 
to iig-iii - Since the upper bound we found for \ \H\\ is a good estimate, this would provide 
us with a best possible condition number. Next, we shall show that one cannot rule out 
a background gauge field configuration for which the lower bound on pprn is saturated. 
These two results put k(H) into a range. The remaining practical question is where in the 
range will k(H) typically be. The answer to this question will depend on details of the 
pure gauge action, and on whether we are dealing with a quenched simulation or with a 
dynamical one. 

Upper bounds for the smallest eigenvalue of H 2 are easily obtained from the variational 
principle. Look at all the block diagonal elements of H 2 : 



Hl, = ^, Hl Xa =c\H 2 w + s 2 s+ 1 ^ Hl, s =H 2 w + s 2 . (50) 



Imagine computing the expectation value of H 2 in a state with an eigenstate of ifjy with 
eigenvalue h in one block, and zero for all other blocks. The variational principle implies 



then Theorem IV. 
Theorem IV: 




(51) 



We assume that h m i n < Then, the minimum is not attained in the ipip block. In the 
other blocks, the best is to take h = h min . Among the XsXs blocks and the (f) s (j) s blocks the 
minimum occurs at s = 1. Thus we find: 

A mi „(# 2 ) < min { ((1 - sin 2 ^) 2 h 2 min + sin 2 |- + ^) , (h 2 min + sin 2 |-) } . (52) 
In the physical relevant cases, /x will be small enough and n will be large enough to give: 

wh - ASr^s- (53) 
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To get a good approximation for the sign function we need n » max \ \ \H\y\ |, tttftjt k so 

I N-"w II J 





• 2 ^ 


ll^ll 2 







\ H w\ 



(54) 



From this we learn that, roughly, the condition number of H cannot be better than the 
condition number of H^y. If we use a low value of n, for which the approximation to the 
sign function is bad, and hence when there are significant violations of chirality beyond 
the explicit mass term fi, the condition number of H might be better than that of Hw- 
This situation however probably defeats the purpose of using the new fermions instead of 
the older, more economical approach of fine tuning the mass term in Hw to small effective 
quark mass. 

On physical grounds it is obvious that there also must exist an upper bound on the 
lowest eigenvalue of if 2 in terms of the "bare" quark mass \l. To see this we recall that the 
basic identity of the entire approach we are focusing on in this paper is 

( H ~^W = ^ 5 + ^e n {H w Y (55) 

The variational principle then gives us Theorem V. 
Theorem V: 



(H 2 ) < A 



1 + n 1 - [i , TT . n2 
■75 H —e n {H w ) 



2 '° 2 

Proof of Theorem V: Let $ = (ipo, 0, 0, ...). On very general grounds, we have 



(56) 



1 



2^ 



> ($|#- 2 |$> > (^ii- 1 !*) 2 = (MH-%^ ). (57) 



Choosing ijj as the eigenstate of (^-^75 + ^-^£n(#wO) of lowest eigenvalue in absolute value 
establishes the inequality ■ 

Note that this inequality is a direct consequence of the existence of some light fermion 
in the theory described by H. Thus, this inequality is very general. 

One cannot exclude backgrounds for which ^(ify^) and 75 have a common eigenvector 
with the eigenvalues of opposite sign. Such eigenvectors should occur in instanton back- 
grounds, for example. For such a background the right hand side of the above equation, 
which is always bounded from below by \i (assume fi > 0, as before) practically saturates the 
bound. In summary, we can say that there are gauge field backgrounds for which we know 
for sure that A mill (if 2 ) < /i 2 . Note that while any eigenvector of H w is also an eigenvector 
of Eooi^Hw) the opposite is not true. 

We see that one cannot hope that the condition number of H would be smaller than either 
the condition number of H\y, or than a number of order -: k(H) > max{K(Hw) , -}■ The 
upper bound on k(H), which we derived previously (equation ([48])), said that the condition 
number of H cannot be worse (larger) than, roughly, K ( Hw > . 

The last question is whether we can imagine a gauge background for which the pessimistic 
(upper) bound (equation (|4"8"D) on k(H) is saturated. The answer is that one can. If the 
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eigenstate of Hw with eigenvalue closest to zero also is an eigenstate of 75 the worst case will 
be realized. In practice however this will happen rarely. The most common low eigenstates 
of H^y, which are the main source of difficulties in simulations, are typically non-degenerate 
and far from chiral (as can be checked by computing the expectation value of 75, which 
would be ±1 for a chiral state). Therefore, in simulations with typical parameters employed 
today, one does not expect to often realize the worst case. 

It is of course important to see how these analytic considerations apply to practical 
simulations. As a toy model we simulated 2-dimensional QED with gauge coupling j3 — 4 
on an 8 x 8 lattice. The main factor governing convergence was the minimal eigenvalue in 
absolute value of H. We found that the upper bound of equation [4^, when viewed as an 
estimate, is overly pessimistic by a factor of order 2. 



VI. DOMAIN WALL FERMIONS 

For domain wall fermions we have, in total, n Dirac fermions labeled by s. The light 
fermion is not as sharply identified as before. Its left handed and right handed components 
reside predominantly at s = 1 and at s = n. The action is ^/D^J, where, adopting notation 
from [HI, we have in four dimensions 



D = D W - P L (M - 1) - P R (M - l)t, P L 



1 + 75 



P 



1-75 



R 



(58) 



Dw is the Wilson Dirac operator with mass m in the range (—2, 0), as before. Dw is unity in 
s-space. M — 1 and (M — l)t have nontrivial action only in s-space, where they approximate 
first order derivatives plus a mass term coupled predominantly to the lightest fermion. 



/ 




M 



10 
10 










1 

0; 



(59) 



D has a hermitian version, i/dwf- 
(S 2 = 1): 



Let the hermitian matrix S produce a flip in s-space 



S 



/o 



1 

\l 



1\ 

1 


0/ 



(60) 



Then, with H 



u = 75-tV, we define 
#dwf = lbSD 



SH W - P L M L - P R M R , 



(61) 



leading to: 
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TJ 2 

-"dwf- 



(62) 



The hermitian mass operators M R = S(M - 1) and M L = S(l - Aft) are related by 
SM R = -M L S. 

Following 0] we can derive a closed formula for det(D^D — z 2 ), but the analysis gets 
complicated and we are not sure that it can be completed. Therefore, we leave this issue 
for the future, and restrict ourselves here to just deriving an upper bound to A min (D^D). 
Since \ m3iX (D^D) would behave similarly to D^Dy/ in one dimension higher, we know that 
\\D\\ will be of the order 10 for QCD [fL4}| . Thus, an upper bound on pprri , will effectively 
provide a best case for n(D^D). Our aim is to show that this best case is similar to what 
we have found for the model we analyzed in the previous sections. 

Let ?/>o be a normalized eigenstate of Hy/ which has the smallest eigenvalue as an eigen- 
state of Construct a trial state in the larger space in which D operates: 



We find: 



H dw ^ = ±Aj£(if£,)tt + 



1 + n 



n 








We now take the norm 



| J D*|| 2 = A min (#^) + 



(l+i^l 2 ±2 I±i^l/2 



n 



n 



This leads to: 



|^|| 2 < KUK) + 



n 



n 



^min(-^iy)- 



Our final bound is given by Theorem VI: 
Theorem VI: 



A m in(-D -D) < 



^min 



i + fi 



i 2 



n 



+ (1 + ^(1-1 



(63) 



(64) 



(65) 



(66) 



(67) 



For domain wall fermions to really represent massless quarks when \i = we need J7J 

» 1. (68) 



n\ 1/2 (H 2 



If we also make the stronger assumption (stronger because we loose nothing by treating only 
the case X m m(Hw) < 1) 
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» 1, 



(69) 



the right hand side of the bound in equation Q87D becomes just \ m i n (Hyy) which is the same 
as in the previous analysis, see equations ([53|,|4]). The difference is in the finite n effects: 
they are larger here. This could be something that works in favor of domain wall fermions. 
Of course, there is no guarantee that the best case, analyzed here, is a good approximation 
to the typical case. One can also analyze trial states with some structure in s-space: 



/ ci-00 \ 
c 2 ^o 



(70) 



V C n 1po } 

The numerical coefficients c s are constrained by 



(71) 



One can optimize the coefficients c s to attain a better bound. 

Let us now argue why one cannot use rigorous methods to get a worst case condition 
number for domain wall fermions that is better than the worst case condition number we 
obtained for overlap fermions implemented by using H (equation fl4*8|)). 

As we learned in the previous section, the key point is that one cannot rule out, for 
arbitrary gauge field backgrounds, the existence of an eigenstate of H w with eigenvalue h 
which also is chiral and for which \h\ is very small. Suppose we have such a state, call it 
ipo- This state is very special: since = Dw, it is a simultaneous eigenstate of Hw and 

Dw- We choose: 



75^o = V>o, H w ijj = D w 4>o = h^ . 



(72) 



Obviously, only a very special gauge background could accommodate a state like this, 
where, in addition, we want \h\ to be very small. Actually, we also want h to be negative. If 
we construct a state ^ out of ipo with coefficients c s as above, we easily see that the action 
of D produces a state of the same structure with only the coefficients c s changed to c' s . The 
action on the coefficients can be immediately read off as: 



/ c i \ 



[1 + h-M] 



V c' n J 



Here, 



/ 




M 



[l + h-M)S 



( a \ 

c 2 

Cn-1 
\ C n J 








M 



( °n \ 

Cn-1 



C 2 

V C! J 



(73) 



-1 1 + h 
\l + h 



-1 l + h\ 
-1 1+h 












(74) 
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The advantage of introducing M is that it is hermitian and obeys 

[1 + h - M] [1 + h - M] f = M 2 (75) 

We need to find the lowest eigenvalue of M for \h\ << 1 and to leading order in /i. We are 
assuming, as usual, that we are close to the chiral limit where \x = 0: 

n\h\ = -nh » 1 (76) 

For /i = we know that the lowest eigenvalue is practically zero [0] (it is responsible for 
| det M\ = \l + h\ n being exponentially small). The eigenstate associated with this eigenvalue 
has the following structure: q = N(l + h) l ~ l (because (1 + h)ci ~ q+i for a very small 
eigenvalue) The normalization N is therefore, at n = oo || 



W = V ££.(i + ft)-"-" = ^ (2 + ft) - (7?) 

Now, first order perturbation theory gives for the smallest (in absolute value) eigenvalue of 
M: 

A min (M) « ~h(2 + h)n, \J \ mm (DW) < |(2 + h)h\n. (78) 

Thus, if a state of type ipQ exists, and under the additional assumption that \h\ « 1, 
we obtain an approximate bound 



< TT^ur ( 79 ) 



\\{DW)-^\\ ~ \\H, 



This leads to essentially the same situation as in the overlap alternative: At the end of 
section V we concluded that one could not hope to be able to prove a better general bound 
on the condition number there than that of equation (pOp. Combining equation ( [T9[ ) with 
the known upper bound | |X?'Z?| | a < (10 + m) jl4| (in four dimensions, with —2 < m < 0), 



produces a bound similar to equation (pEBp. One cannot expect to get a rigorous upper bound 
on the condition number of the Dirac operator employed for domain wall fermions that is 
superior to that obtained for the overlap implemented by H . 



VII. RELATING PARAMETERS 

When comparing the overlap alternative to domain wall fermions one needs some criteria 
to relate the parameters in both models. In a real QCD simulation the criteria should be 
that the parameters be chosen so that similar physics is being described. This is not as 
unambiguous as it first sounds, but clearly beyond this paper. Here we suggest the following 
set of criteria: We need to prescribe relations between the number of extra fields in each 
case, the parameters m in the Wilson Dirac operator in each case and the mass parameters 
\i in each case. 

The easiest is to connect the integers n in the two cases. This we do by requiring the 
error in realizing the sign function in both cases to be the same in the realistic case that Hyy 
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has some small eigenvalues \h\, of order 10 2 or less. The error for domain wall fermions is 
of order 

e -fn. dmf H_ ( 80 ) 
The error in the alternative proposed in this paper is of order 

e -2n|h|_ ( gl ) 

But the number of Dirac fields in the domain wall case is n^mf while in the overlap alternative 
it is n ovp = In + 1 So, for practical values we should take: 

n dm f ~ 2n ovp . (82) 

This ensures similar violations of chirality at \i — for the gauge backgrounds where these 
violations are most significant. This match of the number of extra fields works in favor of 
the overlap alternative. 

To match the parameters m in both cases we focus on the problematic gauge backgrounds 
for which one would need very large numbers of extra fields to reproduce the sign function 
correctly (or else, use projection techniques). We recall that if the parameter m is the same 
in the domain wall and overlap contexts the logarithm of the transfer matrix of domain wall 
fermions and the hermitian Wilson Dirac operator of the overlap both acquire eigenvalues 
very close to zero in the same (bad) gauge configurations jn|. 

What is left is to match the parameters fi. This is more difficult to do in a gauge field 
independent way: Let us adopt the criterion that we want the quark propagator masses to 
be the same in both cases, for the same gauge background, and infinite number of extra 
fields in both cases. The parameter m has already been chosen to be the same. Denoting 
by u the fourth root of the plaquette variable ||16|| , we find, in d dimensions, the following 
"mean-field improved" estimates for the quark masses in each case: 

For domain wall fermions, assuming < /i << 1, we get: 

dwf _ [m + d(l-u)}[2 + m + d(l-u)] 

m phys — {°<JJ 

For the overlap fermions we get ||17|| : 

overlap _ 2/J \m + d(l - u) \ 

m phys ~~ t • l° 4 J 

P y 1 — jl U 

Tree level perturbation theory is obtained by setting u = 1. Above, we assume that —2 < 
m < 0. Setting the two physical masses equal to each other we obtain for small fi 



/^■overlap /^dwf 



m + d(l — u) 

1 + ^ " 



(85) 



In a QCD simulation one may take u w .875, m w —1.8, giving, roughly, /! OV eriap ~ -35/idwf- 
This might indicate an advantage to domain wall fermions, because the slow-down on inver- 
sion of H or D as a result of a small fi might be roughly the same for equal fi parameters. 
In our simulations of the Schwinger model u is closer to unity, u ~ .95 
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VIII. PROJECTION TECHNIQUE 



In practice, most of the numerical problems have to do with the presence of small eigen- 
values to H jy for typical gauge backgrounds. In the older applications of the overlap one used 
a projection method to deal with this difficulty. The projection method trivially extends to 
the overlap alternative presented here. Suppose we have several low states ip a ,H w i(: a = h a ^ a 
with so small \h a \ values that the needed ra's to handle these are too large. Define the 
orthogonal set of projectors 

P a = Ml T, P a='P- (86) 

a 

Redefine Hw in the overlap alternative by: 

H W ^H?=(1- V)H W (1 - V) + $>gn(/gP a . (87) 

a 

This does not change the effective ifi action at n — oo; no new questions about locality 
appear as a result of this replacement. The replacement dramatically reduces the n needed 
to get close to the fi = chiral limit. The shifted states are now perfectly represented 
because for any n we have e n (±l) = ±1. Of course, calculating the P a 's and acting with V 
comes with a cost; it is expected that the cost is bearable because the number of states one 
needed to project out was of the order of 20 in practical simulations carried out so far in 



the quenched approximation using the rational approximation |TT . 

The introduction of projectors complicates the calculation of the force in hybrid Monte 
Carlo, but the complication is manageable |fTS| . Under an infinitesimal change of the back- 
ground gauge fields we have: 

5P a = 5^\ + ^Ml (88) 

We make a phase choice so that P a Sipa = (ip a is normalized; the phase choice elimi- 
nates the component of 8ip a in the direction of ip a and we ignore possible degeneracies for 
simplicity) and then obtain 

<ty« = 5H w ^ a . (89) 

The variation of the projector now follows from that of the state, and the phase choice we 
made no longer has any effect: 

5p a = sH w p a + pjhi . (90) 

The variation of the projectors enters the force always when acting on a vector. To obtain 
the resulting vector we need to compute the action of an inverse on a vector. If the overlap 
Dirac operator is implemented by the rational method of ref J5,ll], the action of c+ 1 Hw 



on the same vector for some constant c has already been evaluated. Using the shift trick 



of reference [0 one can evaluate this new inversion at practically no additional cost in 
operations. 
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Of course, if somebody comes up with a replacement of Hw that has a larger gap around 
zero, it just plugs in simply into the overlap alternative. 

Projection methods when generalized to domain wall fermions need to address the more 
complicated form the transfer matrix Tw- But, if we are willing to depart somewhat from 
the clean structure we have seen above in the overlap case a natural suggestion is to replace 
Hw by H^r in the hermitian domain wall operator i?dwf • We expect this to help because this 
would shift the small energy modes of Hw elsewhere in the spectrum. This is not an exact 
procedure in terms of the transfer matrix T w , but the unit eigenstates of T w (conjugated by 
7 5 ) are also zero eigenstates of H w . Viewing the replacement as a perturbation, we see that 
the troublesome modes are the ones going to be most significantly affected. The main point 
in choosing an approximate projection is to avoid dealing with the nonlocal structure of the 
transfer matrix Tw- Numerical experience from early overlap days teaches us that dealing 
with T w directly is possible JL5|, but it is also rewarding to deal with the sparse matrix H w 



instead. It would be interesting to check this suggestion out in practice. 



IX. ADVANTAGES OF THE OVERLAP ALTERNATIVES 

From the point of view of the projection technique the overlap alternative is certainly 
cleaner. But, if the proposal for domain wall fermions we just made works, we may still feel 
that the domain wall approach and the overlap approach ended up too close to a draw to 
make it worthwhile to investigate the overlap alternatives. 

Let us now turn to what we feel are more clear advantages of the overlap alternatives. 
These advantages have to do with the simplified structure of the effective ip action and 
with the fact that the light fermion is so well identified. Unlike for domain wall fermions, 
where the "wave function" for the light quark penetrates a sizable amount into the extra 
dimension, here the light fermion is fixed as if>. There is no better context to show how this 
would impact numerical QCD than to discuss chiral symmetry at = 0: 

Let us introduce new if)' fields so that the action now involves the non-hermitian overlap 
Dirac operator D : 

V>' = V>, ^ = ^75- (91) 

Introduce yet another fermionic degree of freedom £, but this time an auxiliary one: it has 
no kinetic energy. Integrating out all the fermions but if) and £ leaves the following effective 
action: 

S rf (V/,0 = ^'i±^+fc (92) 

The matrix V n is given by 

V n = ^e n (H w ) (93) 

and obeys: 

V\V n = e 2 n {H w ). (94) 
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At n = oo, for an invertible Hw, we see that V = V^, is unitary and that the ip',^' action 
is given by the overlap Dirac operator, D 

D.= 1 -^. (95) 

It is well known that D^ 1 — 1 is chiral O. 
Define the physical fermion field ip p h by 

V ph = V + £• (96) 

This definition was chosen so that we get, at n = oo, the required subtraction of unity 
(contributed by the £ propagator) from the ip propagator. 

(V^Ph) = D- 1 - 1 = (97) 
Therefore, the effective action for the physical field is 

Sesiipph) = -4> P h 1 ^ ph- (98) 
This effective action is chirally symmetric: The transformation 

V^^Vph, V^ph^Ve^ 75 (99) 
leaves S e e invariant on account of 

1 + V • 1 + V 
e taj 5 _2_ e ^ 5 = _L_ (100) 

A continuum fermion operator with desired chiral properties is now simply transcribed to 
the lattice by replacing the continuum fermionic fields by ^ph- Violations of chirality as a 
result of finite n can be traced quite explicitly. 

Perturbative calculations are an integral part of any procedure that connects numerical 
QCD to physics. In the overlap alternative one can restrict ones attention only to the ip 
fields; there is no need to deal with higher dimensional propagators. One has an explicit, 
relatively simple action, and calculations, although still harder than in the Wilson case, are 
tractable. Indeed, there has been progress on perturbative calculations with overlap fermions 
recently ]2U| (superficially, the treatment of chiral symmetry there may seem somewhat more 
involved then the one presented above, but it is essentially the same as here). As far as we 
know there are no finite n calculations yet, but we expect them to be relatively manageable. 

The parallel calculations in the domain wall case are extremely cumbersome, see for 
example |21|| . In particular, when one deals with questions related to chiral symmetry one 
seems to have to deal with the fermions in the five dimensional "bulk" |22|. This makes 



the finite n effects in the domain wall context harder to estimate quantitatively. All in all, 
when working with domain wall fermions, one pays quite a price for having the two chirality 
components of the light fermion widely separated, and the entire set of extra fermions 
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actively involved in the communication between them, be it because of finite n effects or be 
it because of a topologically nontrivial gauge background. 

Another factor to recall which is in favor of the overlap has to do with the parameter 
m. In the overlap this parameter is theoretically restricted to the range —2 < m < 0. 
When the effects of the noisy gauge background are taken into account one gets effectively 
a positive additive contribution to m (in the mean field approximation it is d(l — u)) which 
forces one towards the —2 end of the range. Still, one does not need to go beyond —2, so 
one is safe even for trivial gauge backgrounds. In the domain wall case one also is driven 
towards —2, but now there is some theoretical worry: For m in the range —2 < m < — 1 
the expression for the transfer matrix Tjy, while still a hermitian matrix, no longer stays 
positive definite for all gauge fields 0. There are gauge backgrounds for which the matrix 
has negative eigenvalues. This raises some concerns about the true phase the lattice model 
is in. Even only the proximity of a phase different from continuum QCD is a source of 
potentially large, undesirable, numerical effects. Theoretically one would like to stay in the 
range — 1 < m < for domain wall fermions, but this cannot be achieved in QCD at present 
typical simulation parameters. In principle, one could get into this range, but the needed 
gauge couplings would have to be impractically weak. Thus, the overlap alternative seems 
safer for the coarse lattices currently employed in numerical QCD. 

Our specific overlap proposal is a very simple implementation of one rational approxima- 
tion to the sign function. Clearly, there are many variations possible and there is room for 
more improvement. In terms of flexibility, the overlap alternative of this paper is superior to 
the domain wall approach; and more flexibility opens more possibilities to increase numerical 
efficiency 

Up to this point our main conclusion is that the overlap approach based on the matrix 
H is superior to the one based on domain wall fermions. Of course, we cannot rule out some 
surprises, so a numerical check in the context of QCD is necessary. 



X. H w OPERATION COUNTS 

Going back to our main motivation, the reason to introduce H in the first place was to 
have a more efficient implementation of the action of the rationally approximated overlap 
Dirac operator. The original direct implementation used a two level nested conjugate gra- 
dient (CG) procedure and a mass-shift trick []5,11,19] which makes operation counts almost 



n- independent. In the inner CG the relevant condition number is that of Hy/ and in the 
outer CG the relevant condition number is roughly -. We see that our worst case bound 
for the condition number relevant to the single step CG for H is roughly the product of the 
previous two condition numbers. Thus, it seems that the new procedure is an order n slower 
than the one it came to improve on. In addition, if one is willing to increase the operation 
count by a factor of 2 in the old procedure, one can eliminate all the storage of the extra 



fields, providing a factor of n saving in memory f23fl ; in a computation limited by memory 
bandwidth rather than processor performance this version might turn out to be the best. 

It is important to mention that these considerations ignore the possibility of precon- 
ditioning the algorithm for inverting H; for example, the structure of H readily admits 
standard red-black preconditioning. Other preconditioning methods might exist, exploiting 
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the rather smooth structure in s-space as evident from the continuum s-limit with no s- 
derivatives. Moreover, the bounds on the condition number are not really saturated very 
often, and super-linear convergence effects in the CG procedure may change the dependence 
of the number of required iterations on the condition number away from the theoretical limit. 
Still, it seems that for very large n the older approach will eventually win. However, both 
the older and newer approach can be improved by projectors, and this will limit the size of n 
one really needs. Also, the mass-shift trick is incompatible with red-black preconditioning, 
so this might work in favor of the new approach. 

To get some feeling for what one would see in practice we again turn to two dimensional 
QED, with a simple plaquette Wilson action at (3 = 4 on an 8 x 8 lattice. We performed 
the calculations necessary to obtain (^ph^ph) using either method and counted the number 
of Hw operations required to reduce the norm of the residual to 10~ 8 . We used n = 20 
but no preconditioning in the ff-algorithm. We did use the mass-shift trick in the older 
algorithm. In neither method did we include projectors. This comparison ended up in favor 
of the older method by a factor of roughly 2.5. This factor could be beaten by red-black 
preconditioning, but we have not tried this out. 

The above direct comparison between the two overlap methods is easy because they 
should produce exactly the same results at the same n and H w , gauge configuration by 
gauge configuration. Essentially, these are just two different algorithms to do the same 
thing. A comparison to domain wall fermions is more difficult, because the differences at 
finite n are more substantial and there is uncertainty about how to match the parameters. 
As far as operation counts go, one action of D in the domain wall case counts roughly as 
ndwf Hw actions. 

A plot of the average number of operations of H\y as a function of mass is shown for 
domain wall fermions, for the direct rational implementation and for the higher dimensional 
implementation of the overlap Dirac operator in Fig. [I]. The data was obtained from a sample 
of 20 configurations. The Wilson mass parameter was set to —1.5 in all cases. Plotted on 
the vertical axis is the number of operations needed for a single inversion of the fermionic 
operator. For the overlap Dirac operator simulations n was set to 20. The plot also includes 
data for domain wall fermions, where n was taken as 40, in accordance with equation (|S2D. 
We used the same parameter fi also for domain wall fermions because the relation in eq. 
(85) holds only for < fi « 1 and becomes totally inadequate at \i = 1. 

XI. CONCLUSIONS 

Let us summarize roughly the situation we were looking at when we began this paper: Al- 
though an approach based on the overlap Dirac operator looked theoretically cleaner, domain 
wall fermions were more attractive numerically. Our analysis has led us to the conclusion 
that there is no evidence that domain wall fermions have even a numerical advantage. 

In all cases we looked at, one faces a problem related to almost zero modes of Hy/. This 
requires large numbers of extra fields in order to preserve chirality. It also affects adversely 
the condition numbers. Whichever method we use, the worst case condition numbers are a 
product of the inverses of two main scale ratios: The first is the scale of the small eigenvalues 
of H^r divided by an upper bound of the order of 5-10 in lattice units. The second scale 



23 



FIG. 1. A comparison of the number of operations of Hw for the inversion of the fermionic op- 
erator in three cases: domain wall fermions, the direct implementation of the rational approximated 
sign function and the higher dimensional implementation of the overlap Dirac operator. 
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ratio is the lattice physical quark mass squared divided by a number of order unity. Each 
small scale ratio slows down inversion independently and the effect compounds in the worst 
case. 

Thus, as far as we can see, at the numerical level, there are no a priori advantages to 
choosing domain wall fermions over overlap fermions in the context of QCD. In both for- 
mulations one faces similar numerical obstacles, and the overlap, to say the least, does not 
fare any worse than domain wall fermions. At the analytical level we are convinced that an 
approach based on the overlap (or any other efficient replacement of the overlap Dirac oper- 
ator that might be found in the future) is superior at presently attainable gauge couplings 
in numerical QCD. Perturbation theory is more transparent to interpret and technically less 
complex in the overlap version. The chirality violating effects associated with the number 
of extra fields are much more explicit and therefore their impact should be easier to trace 
through. 
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